home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dgebal.f < prev    next >
Text File  |  1996-07-19  |  9KB  |  321 lines
  1.       SUBROUTINE DGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
  2. *
  3. *  -- LAPACK routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     September 30, 1994
  7. *
  8. *     .. Scalar Arguments ..
  9.       CHARACTER          JOB
  10.       INTEGER            IHI, ILO, INFO, LDA, N
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   A( LDA, * ), SCALE( * )
  14. *     ..
  15. *
  16. *  Purpose
  17. *  =======
  18. *
  19. *  DGEBAL balances a general real matrix A.  This involves, first,
  20. *  permuting A by a similarity transformation to isolate eigenvalues
  21. *  in the first 1 to ILO-1 and last IHI+1 to N elements on the
  22. *  diagonal; and second, applying a diagonal similarity transformation
  23. *  to rows and columns ILO to IHI to make the rows and columns as
  24. *  close in norm as possible.  Both steps are optional.
  25. *
  26. *  Balancing may reduce the 1-norm of the matrix, and improve the
  27. *  accuracy of the computed eigenvalues and/or eigenvectors.
  28. *
  29. *  Arguments
  30. *  =========
  31. *
  32. *  JOB     (input) CHARACTER*1
  33. *          Specifies the operations to be performed on A:
  34. *          = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0
  35. *                  for i = 1,...,N;
  36. *          = 'P':  permute only;
  37. *          = 'S':  scale only;
  38. *          = 'B':  both permute and scale.
  39. *
  40. *  N       (input) INTEGER
  41. *          The order of the matrix A.  N >= 0.
  42. *
  43. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  44. *          On entry, the input matrix A.
  45. *          On exit,  A is overwritten by the balanced matrix.
  46. *          If JOB = 'N', A is not referenced.
  47. *          See Further Details.
  48. *
  49. *  LDA     (input) INTEGER
  50. *          The leading dimension of the array A.  LDA >= max(1,N).
  51. *
  52. *  ILO     (output) INTEGER
  53. *  IHI     (output) INTEGER
  54. *          ILO and IHI are set to integers such that on exit
  55. *          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
  56. *          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
  57. *
  58. *  SCALE   (output) DOUBLE PRECISION array, dimension (N)
  59. *          Details of the permutations and scaling factors applied to
  60. *          A.  If P(j) is the index of the row and column interchanged
  61. *          with row and column j and D(j) is the scaling factor
  62. *          applied to row and column j, then
  63. *          SCALE(j) = P(j)    for j = 1,...,ILO-1
  64. *                   = D(j)    for j = ILO,...,IHI
  65. *                   = P(j)    for j = IHI+1,...,N.
  66. *          The order in which the interchanges are made is N to IHI+1,
  67. *          then 1 to ILO-1.
  68. *
  69. *  INFO    (output) INTEGER
  70. *          = 0:  successful exit.
  71. *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  72. *
  73. *  Further Details
  74. *  ===============
  75. *
  76. *  The permutations consist of row and column interchanges which put
  77. *  the matrix in the form
  78. *
  79. *             ( T1   X   Y  )
  80. *     P A P = (  0   B   Z  )
  81. *             (  0   0   T2 )
  82. *
  83. *  where T1 and T2 are upper triangular matrices whose eigenvalues lie
  84. *  along the diagonal.  The column indices ILO and IHI mark the starting
  85. *  and ending columns of the submatrix B. Balancing consists of applying
  86. *  a diagonal similarity transformation inv(D) * B * D to make the
  87. *  1-norms of each row of B and its corresponding column nearly equal.
  88. *  The output matrix is
  89. *
  90. *     ( T1     X*D          Y    )
  91. *     (  0  inv(D)*B*D  inv(D)*Z ).
  92. *     (  0      0           T2   )
  93. *
  94. *  Information about the permutations P and the diagonal matrix D is
  95. *  returned in the vector SCALE.
  96. *
  97. *  This subroutine is based on the EISPACK routine BALANC.
  98. *
  99. *  =====================================================================
  100. *
  101. *     .. Parameters ..
  102.       DOUBLE PRECISION   ZERO, ONE
  103.       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  104.       DOUBLE PRECISION   SCLFAC
  105.       PARAMETER          ( SCLFAC = 1.0D+1 )
  106.       DOUBLE PRECISION   FACTOR
  107.       PARAMETER          ( FACTOR = 0.95D+0 )
  108. *     ..
  109. *     .. Local Scalars ..
  110.       LOGICAL            NOCONV
  111.       INTEGER            I, ICA, IEXC, IRA, J, K, L, M
  112.       DOUBLE PRECISION   C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
  113.      $                   SFMIN2
  114. *     ..
  115. *     .. External Functions ..
  116.       LOGICAL            LSAME
  117.       INTEGER            IDAMAX
  118.       DOUBLE PRECISION   DLAMCH
  119.       EXTERNAL           LSAME, IDAMAX, DLAMCH
  120. *     ..
  121. *     .. External Subroutines ..
  122.       EXTERNAL           DSCAL, DSWAP, XERBLA
  123. *     ..
  124. *     .. Intrinsic Functions ..
  125.       INTRINSIC          ABS, MAX, MIN
  126. *     ..
  127. *     .. Executable Statements ..
  128. *
  129. *     Test the input parameters
  130. *
  131.       INFO = 0
  132.       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
  133.      $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
  134.          INFO = -1
  135.       ELSE IF( N.LT.0 ) THEN
  136.          INFO = -2
  137.       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  138.          INFO = -4
  139.       END IF
  140.       IF( INFO.NE.0 ) THEN
  141.          CALL XERBLA( 'DGEBAL', -INFO )
  142.          RETURN
  143.       END IF
  144. *
  145.       K = 1
  146.       L = N
  147. *
  148.       IF( N.EQ.0 )
  149.      $   GO TO 210
  150. *
  151.       IF( LSAME( JOB, 'N' ) ) THEN
  152.          DO 10 I = 1, N
  153.             SCALE( I ) = ONE
  154.    10    CONTINUE
  155.          GO TO 210
  156.       END IF
  157. *
  158.       IF( LSAME( JOB, 'S' ) )
  159.      $   GO TO 120
  160. *
  161. *     Permutation to isolate eigenvalues if possible
  162. *
  163.       GO TO 50
  164. *
  165. *     Row and column exchange.
  166. *
  167.    20 CONTINUE
  168.       SCALE( M ) = J
  169.       IF( J.EQ.M )
  170.      $   GO TO 30
  171. *
  172.       CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
  173.       CALL DSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
  174. *
  175.    30 CONTINUE
  176.       GO TO ( 40, 80 )IEXC
  177. *
  178. *     Search for rows isolating an eigenvalue and push them down.
  179. *
  180.    40 CONTINUE
  181.       IF( L.EQ.1 )
  182.      $   GO TO 210
  183.       L = L - 1
  184. *
  185.    50 CONTINUE
  186.       DO 70 J = L, 1, -1
  187. *
  188.          DO 60 I = 1, L
  189.             IF( I.EQ.J )
  190.      $         GO TO 60
  191.             IF( A( J, I ).NE.ZERO )
  192.      $         GO TO 70
  193.    60    CONTINUE
  194. *
  195.          M = L
  196.          IEXC = 1
  197.          GO TO 20
  198.    70 CONTINUE
  199. *
  200.       GO TO 90
  201. *
  202. *     Search for columns isolating an eigenvalue and push them left.
  203. *
  204.    80 CONTINUE
  205.       K = K + 1
  206. *
  207.    90 CONTINUE
  208.       DO 110 J = K, L
  209. *
  210.          DO 100 I = K, L
  211.             IF( I.EQ.J )
  212.      $         GO TO 100
  213.             IF( A( I, J ).NE.ZERO )
  214.      $         GO TO 110
  215.   100    CONTINUE
  216. *
  217.          M = K
  218.          IEXC = 2
  219.          GO TO 20
  220.   110 CONTINUE
  221. *
  222.   120 CONTINUE
  223.       DO 130 I = K, L
  224.          SCALE( I ) = ONE
  225.   130 CONTINUE
  226. *
  227.       IF( LSAME( JOB, 'P' ) )
  228.      $   GO TO 210
  229. *
  230. *     Balance the submatrix in rows K to L.
  231. *
  232. *     Iterative loop for norm reduction
  233. *
  234.       SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
  235.       SFMAX1 = ONE / SFMIN1
  236.       SFMIN2 = SFMIN1*SCLFAC
  237.       SFMAX2 = ONE / SFMIN2
  238.   140 CONTINUE
  239.       NOCONV = .FALSE.
  240. *
  241.       DO 200 I = K, L
  242.          C = ZERO
  243.          R = ZERO
  244. *
  245.          DO 150 J = K, L
  246.             IF( J.EQ.I )
  247.      $         GO TO 150
  248.             C = C + ABS( A( J, I ) )
  249.             R = R + ABS( A( I, J ) )
  250.   150    CONTINUE
  251.          ICA = IDAMAX( L, A( 1, I ), 1 )
  252.          CA = ABS( A( ICA, I ) )
  253.          IRA = IDAMAX( N-K+1, A( I, K ), LDA )
  254.          RA = ABS( A( I, IRA+K-1 ) )
  255. *
  256. *        Guard against zero C or R due to underflow.
  257. *
  258.          IF( C.EQ.ZERO .OR. R.EQ.ZERO )
  259.      $      GO TO 200
  260.          G = R / SCLFAC
  261.          F = ONE
  262.          S = C + R
  263.   160    CONTINUE
  264.          IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
  265.      $       MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
  266.          F = F*SCLFAC
  267.          C = C*SCLFAC
  268.          CA = CA*SCLFAC
  269.          R = R / SCLFAC
  270.          G = G / SCLFAC
  271.          RA = RA / SCLFAC
  272.          GO TO 160
  273. *
  274.   170    CONTINUE
  275.          G = C / SCLFAC
  276.   180    CONTINUE
  277.          IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
  278.      $       MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
  279.          F = F / SCLFAC
  280.          C = C / SCLFAC
  281.          G = G / SCLFAC
  282.          CA = CA / SCLFAC
  283.          R = R*SCLFAC
  284.          RA = RA*SCLFAC
  285.          GO TO 180
  286. *
  287. *        Now balance.
  288. *
  289.   190    CONTINUE
  290.          IF( ( C+R ).GE.FACTOR*S )
  291.      $      GO TO 200
  292.          IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
  293.             IF( F*SCALE( I ).LE.SFMIN1 )
  294.      $         GO TO 200
  295.          END IF
  296.          IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
  297.             IF( SCALE( I ).GE.SFMAX1 / F )
  298.      $         GO TO 200
  299.          END IF
  300.          G = ONE / F
  301.          SCALE( I ) = SCALE( I )*F
  302.          NOCONV = .TRUE.
  303. *
  304.          CALL DSCAL( N-K+1, G, A( I, K ), LDA )
  305.          CALL DSCAL( L, F, A( 1, I ), 1 )
  306. *
  307.   200 CONTINUE
  308. *
  309.       IF( NOCONV )
  310.      $   GO TO 140
  311. *
  312.   210 CONTINUE
  313.       ILO = K
  314.       IHI = L
  315. *
  316.       RETURN
  317. *
  318. *     End of DGEBAL
  319. *
  320.       END
  321.